######### This script will remove small patches from a species realized distribution ASCII
#########

#### Establish base directory

base.dir = '/home1/99/jc152199/MAXENT/'
setwd(base.dir)

#### Load library

library('SDMTools')

##### Need to loop through ASCIIs for multiple species
##### First identify ASCII directory

ascii.dir = '/home1/99/jc152199/MAXENT/output/'

#### Directory to write files to

out.dir = '/home1/99/jc152199/skinkpatchstatrealizeddist/'

##### Identify species

spp = dir(ascii.dir)

#### Remove directories which aren't skinks

spp = spp[c(1,13:15,17:19)]

###### Now loop through and identify individual ASCII for each species
##### Don't forget, this needs to be done for both the microCLIM, expCLIM, and BIOCLIM distributions

for (s in spp)

	{
	
	#### Species loop
	
	for (m in c('microCLIM','BIOCLIM'))
	
		{
		
		#### Model loop
		
		#### Read in an ASCII of realized distribution
		
		tasc = read.asc.gz(paste(ascii.dir,s,'/',m,'/output/',s,'_',m,'_Realized.asc.gz',sep=''))
		
		#### Keep an unchanged copy of tasc called base.asc
		
		base.asc =  tasc
		
		#### Convert to a binary matrix
		#### Don't need to apply MaxEnt thresh because that was done before producing the realized distribution
		
		tasc[which(tasc>0)]=1
		
		#### Run class-stat on tasc to determine total area for the species distribution
		
		cstat = ClassStat(tasc,cellsize=250,bkgd=NA,latlon=TRUE)[2,]
		
		### Extract the total area of the species realized distribution
		
		s.area = cstat[1,3]
		
		#### Run ConnCompLabel to identify individual patches
		
		pasc = ConnCompLabel(tasc)
		
		#### Run Patch Stat on pasc
		
		pstat = PatchStat(pasc, cellsize=250, latlon=TRUE)
		
		### Write out pstat to the out.dir
		
		write.csv(pstat, file=paste(out.dir,s,'_',m,'_PatchStatSummary.csv',sep=''))
		
		}
		
	cat('\n',s,' - Completed\n',sep='')	
		
	}
	
#### Done

### Now perform t-tests and compare variances between models within species

outdata = NULL

for (s in spp)

	{
	
	### Read in patch stat summaries for a species
	
	acps = read.csv(paste(out.dir,s,'_microCLIM_PatchStatSummary.csv',sep=''), header=T)
	bcps = read.csv(paste(out.dir,s,'_BIOCLIM_PatchStatSummary.csv',sep=''), header=T)
	
	#### Remove first row (stats for no occurrences) and first column (row names)
	
	acps = acps[-1,-1]
	bcps = bcps[-1,-1]
	
	#### Determine degrees of freedom for both datasets (n-1)
	
	acdof = nrow(acps)-1
	bcdof = nrow(bcps)-1

	#### Now begin looping through individual stats
	
	for (stat in names(acps)[-1])
	
		{
		
		### Select a single statistic to work with
		
		acss = acps[,which(names(acps)==stat)]
		bcss = bcps[,which(names(bcps)==stat)]
		
		### Now calculate the variance of the two datasets
	
		acvar = var(acss)
		bcvar = var(bcss)
		
		#### Calculate the F-Ratio (larger variance value is always the numerator)
		
		if (acvar>bcvar)
		
			{
			
			#### Calculate the Critical F-Value
			#### This means the variance of one data set must be CritF times greater than the other before we can determine the variances are significantly different
			#### Homogeneous variance is an important assumption for the t-test
	
			CritF = qf(.975,acdof,bcdof)
		
			### Calculate the F Statistic
		
			FRatio = acvar/bcvar
			
			### Calculate the p-value for an F Statistic with x degrees of freedom
		
			Fp = 2*(1-pf(FRatio,acdof,bcdof))
			
			}
			
		if(bcvar>acvar)	
			
			{
			
			#### Calculate the Critical F-Value
			#### This means the variance of one data set must be CritF times greater than the other before we can determine the variances are significantly different
			#### Homogeneous variance is an important assumption for the t-test
	
			CritF = qf(.975,bcdof,acdof)
			
			### Calculate the F Statistic
			
			FRatio = bcvar/acvar
			
			### Calculate the p-value for an F Statistic with x degrees of freedom
		
			Fp = 2*(1-pf(FRatio,bcdof,acdof))
			
			}
		
		### If variances are significantly different, complete t-test and bind data
		
		if(FRatio<abs(CritF))
		
			{
			
			### t-test
			
			tt = t.test(acss,bcss)
			
			if(tt$p.value<=.05)
				
				{
			
				t.out = data.frame(spp=s, stat = stat, equal.variance = 'YES', CritF=round(CritF,2), FCalc=round(FRatio,2), Fp = round(Fp,2), test='t-test',sig.diff='YES', Wp=NA, tp=round(tt$p.value,2))
		
				}
				
			if(tt$p.value>.05)
			
				{
				
				t.out = data.frame(spp=s, stat = stat, equal.variance = 'YES', CritF=round(CritF,2), FCalc=round(FRatio,2), Fp = round(Fp,2), test='t-test',sig.diff='NO', Wp=NA, tp=round(tt$p.value,2))
		
				}
			
			}		
		
		### If variances are significantly different, complete a Wilcox test then bind data
		
		if(FRatio>=abs(CritF))
		
			{
			
			#### Perform wilctest
			
			wilc.test = wilcox.test(acss,bcss)
			
			if(wilc.test$p.value<=.05)
			
				{
		
				t.out = data.frame(spp=s, stat = stat, equal.variance = 'NO', CritF=round(CritF,2), Fp = round(Fp,2), FCalc=round(FRatio,2), test='Wilcox',sig.diff='YES', Wp=wilc.test$p.value, tp=NA)
		
				}
				
			if(wilc.test$p.value>.05)

				{
				
				t.out = data.frame(spp=s, stat = stat, equal.variance = 'NO', CritF=round(CritF,2), Fp = round(Fp,2), FCalc=round(FRatio,2), test='Wilcox',sig.diff='NO', Wp=wilc.test$p.value, tp=NA)
				
				}
		
			}
		
		### Bind data
		
		outdata = rbind(outdata,t.out)
		
		#### Report Progress
		
		cat('\n',s,' ',stat,' Completed\n',sep='')
		
		}
		
	#### Report progress

	cat('\n',s,' All Stats Completed\n',sep='')
		
	}
	
#### Close loop

#### Still not too certain about the t-test, specifically the calculation of the F-Ratio, and where which mean/var value goes where in the formula

### Write outdata

write.csv(outdata, file=paste(out.dir,'SigTestForAllSkinksPatchStat.csv',sep=''), row.names=F)

### Done

nonequalvar = (length(which(outdata$equal.variance=='NO'))/nrow(outdata))*100

nonsig = (length(which(outdata$sig.diff=='NO'))/nrow(outdata))*100

sig = (length(which(outdata$sig.diff=='YES'))/nrow(outdata))*100

#################################################################
#################################################################

### Aggregate by species and stat type and count the significant differences

sigbyspp= NULL

for (s in spp)

	{
	
	tdata = outdata[which(outdata$spp==s),]
	
	tsig = tdata[which(tdata$sig.diff=='YES'),]
	
	t.out = data.frame(spp=s, n.sig = nrow(tsig))
	
	sigbyspp = rbind(t.out,sigbyspp)
	
	}
	
#### Done

sigbystat = NULL

for (st in unique(outdata$stat))

	{
	
	tdata = outdata[which(outdata$stat==stat),]
	
	tsig = tdata[which(tdata$sig.diff=='YES'),]
	
	t.out = data.frame(stat=st,n.sig = nrow(tsig))
	
	sigbystat = rbind(t.out, sigbystat)
	
	}
	
### Done


	
	